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The Bedraiga's skink gecko, Teratoscincus bedriagai, is one of the fairly known lizards of 
Iran. Phylogenetic relationships among five distant populations of ‘Teratoscincus 
bedriagai were estimated based on the mitochondrial COI partial fragment. This study 
highlights very low genetic divergences among examined populations despite the 
noticeable geographic distance between them. We assumed that the current genetic 
structure of haplotypes is shaped as a consequence of climate changes and glaciations in 
the Quaternary. 
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INTRODUCTION 
The frog-eyed geckos of the genus Teratoscincus Strauch, 1863, have a wide distribution range along 
the arid belt of Middle and Central Asia (Szczerbak & Golubev, 1996; Anderson, 1999; Macey ef a/, 
1999, 2005). Iran with five out of eight described species (Uetz & HoSek, 2017), reflects a unique 
biogeographic condition among the inhabited countries. However, for quite a while, the Iranian 
populations of the genus Teratoscincus haven't been subject to taxonomic research, with the recent 
exception resulting in the description of new taxa (Nazarov ef a/, 2017; Akbarpour ef a/., 2017). 

Teratoscincus bedriagai is one of the slightest well-known reptiles that occurred in Iran 
(Szezerbak & Golubev, 1996; Anderson, 1999). There have been few studies on T\bedriagai in Iran 
over the previous few years. Only morphological data (Hojati e¢ a/, 2009) and something on its 
biology (Hojati e¢ a/, 2014), both from a confined range in Semnan Province are available. 
Therefore, there is yet no sufficient data on the geographic variation in morphological characters in 
T. bedriagai, and also its genetic structure in Iran. One of the most recent molecular studies show that 
the phylogenetic position of T. bedriagai is poorly resolved among other species of Teratoscincus 
(Nazarov ef a/., 2017). 

This study provides a preliminary assessment of the intraspecific variability among distant 
populations of T. bedriagai in Iran. For this purpose, phylogenetic analyses were conducted based on 
the partial COI sequences from five populations within its distribution range in Iran. 
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FIGURE 1. Map of the sampling localities for Teratoscincus bedriagai in Iran. Black lozenges represent 
the sites where specimens of present study coming from while black triangle represent a locality near 
Gonabad, Razavi Khorasan Province, used in Nazarov ef a/ (2017). For more details of numbers, 
refer to Table 1. 


MATERIAL AND METHODS 
As a source of DNA, we used 15 ethanol-preserved tissue samples from five geographically distant 
populations of Teratoscincus bedriagai in Iran (Table 1 & Fig. 1). Voucher specimens were deposited in 
the Herpetological collections of Hakim Sabzevari University, Sabzevar (ERP) and Shahrekord 
University, Shahrekord (HAC). 

We performed DNA extraction using Aron-Gene Kit (Aron-Gene, Isfahan, Iran) following 
the instructions provided by the manufacturer. A fragment of the mitochondrial cytochrome C 
oxidase subunit I (COI) was amplified using combinations of primers L 1498 (Folmer et a/, 1994), 
RepCOI F and RepCOI R (Nagy éf a/, 2012). Standard PCR reaction was done as described by 
Khosravani e¢ a/. (2017). The amplified DNA was then visualized and assessed on a 1.5% agarose 
gel. Positive PCR products were sent for sequencing to Macrogen Inc. (Seoul, South Korea, 
http://www.macrogen.com). All newly obtained sequences were submitted in GenBank under 
accession numbers MH476256 - MH476270 (Table 1). Published sequences of Teratoscincus bedriagai 
(MF573787- MF573790) were downloaded from NCBI (http://www.ncbi.nlm.nih.gov) and 
included in the final dataset. For outgroups, we used T. Reyserlingii (MF573791 and MF573792) and 
T. microlepis (MF573800 and MF573798) which are closely related to T. bedriagai (Nazatov ef al, 
2017). 
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TABLE 1. Information on the Teratoscincus bedriagai specimens used in the study, with locality and 
Genbank accession numbers. 


Voucher = Haplotype : ; Genbank 
NO specimens Locality Neaier Latitude /longitude Niue 
Semnan Province, Kaal-e H1 
ERP 5470 Shour/60 Kirt se Bianowand 35.88, 56.33 MH476256 
Semnan Province, Kaal-e H2 
1 ERP 5471 Shour-60 Knit B.of Bianomatd 35.88, 56.33 MH476257 
Semnan Province, Kaal-e H3 
ERP 5472 Shout 60 Km Hor Bianomand 35.88, 56.33 MH476258 
Razavi Khorasan province, on H9 
ERP 3868 _ the road of Niaz-Abad to 34.1818, 60.31088 MH476264 
Moussa-Abad 
2 Razavi Khorasan province, on H10 
ERP 3869 the road of Niaz-Abad to 34.1818, 60.31088 MH476265 
Moussa-Abad 
Razavi Khorasan province, H14 
ERP 6042 reancdiNe anne --- MH476270 
South Khorasan Province- On H11 
ERP 5974 _ the road of Nehbandan to Deh- 31.4582, 59.6094 MH476266 
Salm- 45 km E of Deh-Salm 
South Khorasan Province- On H12 
3 ERP 5975 the road of Nehbandan to Deh- 31.4582, 59.6094 MH476268 
Salm- 45 km E of Deh-Salm 
South Khorasan Province- On H11 
ERP 5977 the road of Nehbandan to Deh- 31.4582, 59.6094 MH476267 
Salm- 45 km E of Deh-Salm 
Sistan and Baluchistan H4 
HAC 172 Province, NW Zahak, near the Shido Biee MH476259 
ca 61.631615 
Khamak village 
Sistan and Baluchistan H5 
HAC 173 Province, NW Zahak, near the pede MH476260 
aes 61.631615 
Khamak village 
Sistan and Baluchistan H6 
4 HAC 174 Province, NW Zahak, near the es MH476261 
; 61.631615 
Khamak village 
Sistan and Baluchistan H7 
HAC 175 Province, NW Zahak, near the Hagens MH476262 
; 61.631615 
Khamak village 
Sistan and Baluchistan H8 
HAC 176 Province, NW Zahak, near the Bs MH476263 
. 61.631615 
Khamak village 
5 ERP 831 Indefinite locality from Razavi H13 —_ MH476269 


Khorasan province 


Genetic distances within and between populations were calculated using Mega 7.0 (Kumar ef a/, 
2016). The Median-joining algorithm (Bandelt e¢ a/, 1999) implemented in Network 5.0.0.3 
(www.fluxus-engineering.com) was used to depict divergence among haplotypes. The DNA 
sequences were divided initially into three partitions, the first, second and third codon positions of 
COI. A greedy search with PartitionFinder2 (Lanfear et a/., 2016) was conducted to select optimal 
partitioning schemes and the best-fit evolutionary models based on the Akaike information criterion 
(AIC) allowing for separate estimation of branch lengths for each partition. The suggestion of 
PartitionFinder2 analysis is using the predefined partitions except for the second codon positions 
which should be combined in a single partition with first codon positions. Bayesian inference (BI) of 
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phylogeny was performed using MrBayes v.3.2.5 (Ronquist ef a/, 2012). Two simultaneous and 
independent runs consisting of eight MCMC chains in MrBayes were run for 10,000,000 generations 
with default priors, trees sampled every 1,000 generations and separate estimation of parameters for 
individual partitions under default heating using best-fit models as suggested by PartitionFinder2 
(HKY-+I for first and second codon positions as one subset and GTR for third codon positions). 
We assessed convergence of the runs by discarding the 25% of the initial samples as burn-in and a 
majority rule consensus tree was generated from the remaining post-burn-in trees. Maximum 
Likelihood (ML) analysis was performed using RaxML GUI v. 0.95 (Silvestro & Michalak, 2012) 
with ‘GTR’ option and 1000 bootstrap replicates (partition set by codon and 10 independent runs). 


RESULTS 

The dataset used for the phylogenetic analyses consisted of 608 nucleotides, there were 435 
invariable sites and 173 polymorphic sites, of which 141 sites were phylogenetically informative and 
32 were singletons. Our dataset could be translated into amino acid sequences with the vertebrate 
mtDNA genetic code. Nucleotide composition analyses revealed 17 haplotypes, haplotype diversity 
of 0.9883 and the average of 0.496 for GC content for Teratoscincus bedriagai. 

Uncorrected genetic p-distances between five populations of T. bedriagai and outgroup taxa 
are shown in Table 2. The minimum genetic distance between pairs of populations referred to 
Gonabad and Biarjomand was 0.009 while the maximum divergence between Deh-Salm and Niaz- 
Abad was 0.024 (Table 2 and Fig.1). The genetic distance between T. bedriagai with T. microlepis and 
T. Reyserlingi was almost the same (0.171 and 0.166 respectively). The tree inferred from the Bayesian 
analysis was similar to those constructed using Maximum Likelihood (not shown) generally (Fig. 2). 
In both trees, IT. bedriagai recovered as a monophyletic group relative to the outgroup taxa, T. 
microlepis and T. Reyserlingit, with high statistical supports, PP (1) and ML (100%) (Fig. 2). Included 
populations of T. bedriagai show relatively high genetic homogeneity. Nevertheless, a clade including 
individuals from Sistan, Deh-Salm and one sample from Niaz-Abad (Fig. 1) was discovered with 
strong support (100%) in ML tree, albeit not find in Bayesian framework. Some degrees of 
polytomous relationships among T. bedriagai individuals revealed in both analyses. Median-joining 
haplotype network of 17 haplotypes is illustrated in Figure 3. 


TABLE 2. Uncorrected genetic p-distances among five populations Teratoscincus bedriagai and 
outgroup taxa. Since one individual (ERP 831) comes from an indefinite locality in Razavi Khorasan 
province, it was not included. 


1 2 3 4 5 6 

1 T. Reyserlingit 

2 T. microlepis 0.146 

3 T. bedriagai (Biarjomand) 0.167 0.173 

4 T. bedriagai (Niaz-Abad) 0.171 0.174 0.019 

5 T. bedriagai (Deh-Salm) 0.167 0.173 0.021 0.024 

oF bedriagai (Khamak) 0.166 0.173 0.017 0.021 0.015 

7 _T. bedriagai_ (Gonabad) 0.160 0.166 0.009 0.016 0.015 (0.012 
DISSCUSSION 


The genetic variability between five populations of Teratoscincus bedriagai was analyzed using mtDNA 
COI partial sequences. Although T. bedriagai as a monophyletic group was recovered in our analysis. 
However, our data added nothing more to the previous finding on its phylogenetic position among 
other species of the Genus Teratoscincus (Nazarov et al, 2017). Inferring phylogenetic relationships 
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based on gene trees is one of the most reason for current obstacle, which most likely might be 
improved and shed new light on the issue by taking advantages of multilocus phylogeny (e.g. 
Martens ef a/., 2008; Liu e¢ a/., 2017). 

Our findings are compatible with the current taxonomy of T. bedriagai as a valid species with 
numerous populations in Iran (Uetz & HoSek, 2017). This conclusion is derived from little genetic 
differences between T. bedriagai populations in Iran (Table 2), together with their phylogenetic 
relationships with other species of the Genus Tera/oscincus (Fig. 2).According to Star-like topology of 
Median-joining network (Fig.3) and Gonabad's geographic position (H15 likely ancestral haplotype) 
located in the center of the rest of populations, a recent demographic expansion could be a plausible 
scenario. However, the spatial distribution of the haplotypes (Fig. 3) could be the consequences of 
multiple colonization and extinction events recently. It seems that current distribution pattern of T. 
bedriagai is shaped by a huge influx of climate changes and glaciations in the Quaternary (Taberlet et 
al., 1998; Hewitt, 2004; Popov e¢ a/, 2006), but additional analyses of more populations of T. bedriagai 
are necessary to testify this hypothesis. 


T. microlepis 


TF. microlepis 


T. keyserlingii 


T. keyserlingii 


FIGURE 2. The Phylogenetic tree resulting from a Bayesian framework reconstruction based on 608 
base pairs of the COI sequence. Bootstrap support values (left) from maximum likelihood analysis 
and posterior probabilities (right) from Bayesian inference are indicated at the nodes. The tree was 
rooted by Teratoscincus Reyserlingii. 
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FIGURE 3. A Median-joining network connecting 17 haplotypes found in the COI sequence data set. 
Median vector shown with orange nodes, the branches length is indicated the numbers of 
mutational steps joining the haplotypes and the size of each circle is proportional to the number of 
haplotypes represented. Biarjomand (H1: ERP 5470, H2: ERP 5471, H3: ERP 5472); Khamak 
village (H4: HAC 172, H5: HAC 173, H6: HAC 174, H7: HAC 175, H8: HAC 176); Niaz-Abad (H9: 
ERP 3868, H10: ERP 3869, H14: ERP 6042); Deh-Salm (H11: ERP 5974 and ERP 5977, H12: ERP 
5975); Gonabad (H15: 1913 and 1914, H16: 039a, H17: 1912) and H13: ERP 831. For more 
information, refer to Table 1. 
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